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1.  Project  Summary 

Project  Goal 

The  two  primary  objectives  of  the  “Foundations  for  MEMS  Synthesis”  project  were  to 
develop: 

•  A  VLSI-like  Hierarchy  for  MEMS 

•  MEMS  Synthesis  Flow 


Our  technical  approach  comprised  four  major  tasks  that  focus  on  the  circuit-level  design  of 
MEMS.  The  tasks  are  partitioned  by  placement  in  the  synthesis  flow:  structural  synthesis,  circuit/ 
layout  synthesis,  shape  and  process  optimization,  hierarchical  representation  and  performance  veri¬ 
fication,  and  technology  integration  and  transfer.  Figure  1  illustrates  the  relationship  between  these 
research  areas.  Each  of  the  participating  investigators  in  the  tasks  developed  prototypical  implemen¬ 
tations  of  point  design  tools  that  fit  within  the  synthesis-centric  MEMS  design  methodology.  Across 
all  the  tasks,  these  implementations  support  the  overall  design  flow  through  use  of  a  common  hier¬ 
archical  representation  of  MEMS  and  use  of  a  common  set  of  components. 

When  the  project  began,  the  primary  design  representations  for  MEMS  were  layout  and  3D 
solid  modeling  with  customized  meshes  for  each  finite/boundary  element  simulator.  The  choice  of 
the  layout  representation  was  itself  a  challenge. 
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FIGURE  1.  Relationship  between  project  tasks  for  the  MEMS  synthesis  flow. 
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from  the  commonplace  layout  representation  to  a  schematic  representation.  However,  the 
list  of  elements  in  the  component  database  remained  the  same.  This  project  has  therefore  partici¬ 
pated  in  two  revolutions  in  design  representations  for  MEMS  during  the  course  of  the  project. 

Development  in  this  project  of  the  behavioral  simulation  environments  at  Carnegie  Mellon 
and  U.C.  Berkeley  have  enabled  a  revolution  in  the  design  representation.  Specifically,  the  introduc¬ 
tion  of  the  design  representation  abstractions  at  the  schematic  and  layout  level  instead  of  the  3D 
solid  model  or  mesh  level  addresses  the  first  of  the  two  primary  goals  of  this  project.  These  levels  of 
abstractions  now  enable  a  MEMS  designer  to  iteratively  design  a  complete  MEMS  sensor  or  actua¬ 
tor  system,  instead  of  just  a  focused  portion  of  the  transducer  element. 

The  second  goal  of  the  project,  a  synthesis-based  design  flow,  allows  expert  MEMS  design¬ 
ers  to  capture  and  encode  design  constraints  and  enables  novice  MEMS  designers  to  develop 
MEMS  designs  that  have  a  higher  chance  of  working  during  the  first  pass  at  fabrication.  During  the 
project  we  demonstrated  both  the  individual  tools  in  this  flow,  as  well  as  some  interactions  between 
the  tools. 

The  impact  of  these  new  levels  of  abstraction  and  the  new  design  methodology  is  the  ability 
of  the  designer  to  directly  link  the  impact  of  MEMS  device  level  innovations  with  quality  improve¬ 
ments  in  the  final  system  design  as  well  as  enable  the  top-down  design  of  MEMS 
[1][13][14][25][29][33]. 

2.  Capabilities  Developed  During  the  Project 

The  set  of  capabilities  developed  during  the  project  to  support  the  hierarchical  design  repre¬ 
sentation  and  synthesis  flow  are  detailed  for  each  task.  A  tutorial  summary  of  much  of  this  work 
will  be  published  in  IEEE  Transactions  on  CAD  [38].  A  draft  preprint  of  this  paper  is  attached  to  the 
report. 

Task  1:  Structural  Synthesis  (CMU) 

Structural  synthesis  involves  choosing  a  design  topology  that  meets  the  required  engineering 
specifications  for  the  MEMS  device.  In  this  project  we  developed  a  shape-grammar  based  topology 
representation  for  MEMS  [20][22][30],  an  A-Design  based  search  algorithm  to  choose  amongst 
competing  topologies  [16][21],  and  a  hierarchical  method  for  rapid  evaluation  of  the  performance 
of  candidate  design  topologies.  The  resulting  integration  of  these  three  concepts  were  used  to  dem¬ 
onstrate  a  methodology  to  configure  new  resonator  and  accelerometer  structures  automatically. 
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FIGURE  2.  Structural 
Synthesis  using  A-Design 
proposes  MEMS 
topologies  that  transform 
the  user  specified  input 
functional  domain  into  the 
user  specified  output 
functional  domain. 
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By  representing  the  desired  input  and  output  energy  domains,  and  using  a  component  data¬ 
base  that  provides  elements  that  can  translate  between  various  energy  domains,  the  A-Design  based 
search  algorithm  identifies  candidate  topologies  that  meet  the  designers  requirements.  The  use  of  a 
shape-grammar  for  topology  representation  leads  to  the  encoding  of  the  MEMS  layout  (i.e.,  shape), 
allowing  the  search  algorithm  to  directly  interface  with  the  common  hierarchical  representation  for 
the  synthesis-centric  design  flow.  The  hierarchical  evaluation  method  is  based  on  behavioral  simula¬ 
tion  environments  detailed  in  Task  4,  and  allows  the  A-design  search  engine  to  have  a  range  of  per¬ 
formance  evaluators  ranging  from  fast  with  less  accuracy  (for  the  early  stages  of  A-design)  to  slow 
and  accurate  (for  the  final  iterations  in  A-design).  The  example  in  Figure  2  shows  acceleration  input 
and  voltage  output  respectively  and  shows  that  the  structural  synthesis  methodology  led  to  topolo¬ 
gies  that  are  different  from  conventional  manual  design. 

Task  2:  Circuit/Layout  Synthesis  (CMU) 

The  circuit/layout  synthesis  algorithms  developed  during  this  project  can  be  employed  to 
optimally  size  a  fixed  topology  to  meet  the  user  specified  engineering  performance  parameters.  The 
fixed  topology  could  be  generated  by  the  automated  methodology  in  Task  1  or  through  manual 
design.  Three  topology-specific  synthesis  modules  were  developed  in  this  task:  a  folded-flexure  res¬ 
onator  in  the  MUMPS  process  [2][11][23][37],  a  lateral  accelerometer  in  the  MUMPS  process 
[26]  [27],  as  well  as  an  accelerometer  in  the  CMOS-MEMS  process  [36].  The  CMOS-MEMS  accel¬ 
erometer  ties  into  the  joint  DARPA  MTO  MEMS  and  Composite  CAD  funded  project  on  “Inte¬ 
grated  MEMS  Inertial  Measurement  Unit.”  These  synthesis  modules  are  based  on  lumped 
parameter  models  that  map  the  design  parameters  defining  the  layout  geometry  into  the  design  per- 
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FIGURE  3.  Synthesized  layouts  of  (a)  family  of  resonators  optimized  for  area,  voltage  and  thickness; 
MUMPS  accelerometers  under  (b)  open-loop  and  (c)  closed-loop  operation  optimized  range;  and  (e) 
CMOS-MEMS  accelerometers  compared  to  (f)  manual  design. 


formance  space,  and  generate  an  optimization-based  layout.  These  models  incorporate  manufactur¬ 
ing  variations,  such  as  those  affecting  cross-axis  sensitivity  in  accelerometers,  as  well  as  the 
capacitive  load  and  noise  parameters  of  the  circuits  that  condition  the  signal  output  from  the  MEMS 
sensor.  Example  synthesized  layouts  from  each  of  the  modules  are  shown  in  Figure  3. 

We  have  demonstrated  that  such  synthesis  modules  enable  designers  of  integrated  MEMS  to 
synthesize  layouts  that  meet  their  design  needs  on  demand.  These  modules  demonstrate  the  crucial 
MEMS  design  trade-offs  for  a  specific  application  to  novice  MEMS  designers.  Cronos  Integrated 
Microsystems  is  using  a  web-based  interface  to  these  synthesis  tools  in  their  short  course  on 
MEMS.  Microcosm  Technologies,  Inc.  is  licensing  the  technology  developed  during  this  task  with 
plans  for  offering  it  as  a  teaching  aid.  The  topology  specificity  of  the  synthesis  modules  and  the  cur¬ 
rent  lack  of  automation  in  macromodel  generation  prevents  a  more  general  tool  at  this  time. 

To  overcome  this  limitation,  progress  has  been  made  to  integrate  the  optimization  frame¬ 
work  developed  in  Task  2  with  the  behavioral  simulation  methodology  in  Task  4.  This  integration 
will  allow  any  user  comfortable  with  the  simulation  methodology  to  use  the  circuit/layout  synthesis 
technologies  for  automated  design.  Additionally,  we  fabricated  four  optimized  CMOS-MEMS 
accelerometers  to  experimentally  verify  the  models  in  the  synthesis  module  (previous  experimental 
fabrication  of  MUMPS  resonators  showed  that  the  models  compared  to  within  10%  of  the  fabri¬ 
cated  devices).  Testing  of  these  devices  will  occur  under  another  DARPA-MTO  project. 
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Task  3:  Shape  and  Process  Optimizations  (MIT,  UPenn) 

The  MIT  efforts  focused  on  developing  fast  algorithms  for  computing  sensitivities  in  numer¬ 
ical  simulation,  to  provide  search  directions  for  the  shape  optimization  efforts  at  UPenn,  as  well  as 
for  developing  a  framework  for  process  optimization.  Fast  computations  for  electrostatic  sensitivi¬ 
ties  were  demonstrated  in  the  first  year  [6],  and  a  rigid  elastic  formulation  coupled  with  multi-level 
Newton  solution  subsequently  demonstrated  additional  speedups  [9].  Optimization  algorithms  for 
obtaining  the  shape  of  compliant  structures  to  meet  the  desired  force-deflection  characteristics  were 
demonstrated  at  UPenn  [17] [19].  As  an  extension  of  the  static  force-deflection  problem,  the  struc¬ 
tural  synthesis  for  dynamic  specifications  such  as  the  (eigen)  modeshapes  of  a  ring  gyroscope  were 
demonstrated  [18].  The  third  subtask  on  process  optimization  frameworks,  at  MIT,  led  to  the  devel¬ 
opment  of  process  representation  and  repository  methods  for  capturing  manufacturing  variation 
(implemented  using  the  “Semiconductor  Process  Representation”  standard);  propagation  of  process 
variation  to  structure  variation;  and  analysis  of  process  variation  on  device  performance.  This  task 
was  linked  with  the  rest  of  the  project  in  several  ways:  one  of  the  static  force-deflection  examples  at 
UPenn  involved  the  maximizing  the  mechanical  displacement  of  an  inertial  mass  using  a  compliant 
mechanism  to  link  with  the  accelerometer  synthesis  efforts  in  Task  2.  Also,  the  resonator  synthesis 
module  of  Task  2  was  modified  to  link  with  the  process  variation  framework  developed  in  Task  3  to 
trade-off  between  resonator  area  and  frequency  variation. 

Task  4:  Hierarchical  Representation  and  Performance  Verification  (CMU,  UC  Berkeley) 

The  Carnegie  Mellon  and  UC  Berkeley  efforts  independently  demonstrated  the  concept  of 

composing  MEMS  devices  from  a  library  of  parameterized  behavioral  models.  Additionally, 
through  the  gyroscope  canonical  problem,  a  multi-level  hierarchical  design  methodology  was  dem¬ 
onstrated.  The  tools  developed  in  this  sub-task  have  created  a  powerful  new  schematic  design-entry 
mode  for  MEMS  as  an  alternative  to  solid-modeling  and  layout. 

By  developing  an  initial  hierarchical  representation,  and  a  library  catalog  of  “atomic”  ele¬ 
ments  found  in  suspended  MEMS,  and  subsequently  developing  behavioral  models  of  these  ele¬ 
ments,  the  performance  verification  of  resonators  [8],  accelerometers  [10],  vibratory  rate 
gyroscopes  [7],  and  micropositioners  [35]  were  demonstrated.  Subsequent  improvements  in  the 
hierarchical  representation  simplified  the  user  interface  between  the  designer  and  the  simulation 
environment,  and  also  led  to  a  speed-up  in  simulation  times.  Additional  improvements  in  the  ele¬ 
mental  behavioral  models,  and  in  models  of  higher  level  elements  (such  as  the  comb-drive)  led  to 
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FIGURE  4.  (a)  High  level  NODAS  schematic  of  canonical  gyroscope  consisting  of  four  spring  and  comb 
sections  a  plate  mass  and  several  circuits,  with  a  (b)  hierarchical  lower-level  schematic  of  one  of  the  4  spring 
and  comb  sections. 

improved  simulation  accuracy  as  well  as  simulation  performance.  The  current  Cadence-based  Carn¬ 
egie  Mellon  tool  (NODAS)  [8]  [10]  [24]  [28]  [31]  [32]  [34]  and  MATLAB-based  UC  Berkeley  tool 
(SUGAR)  [12] [15]  represent  the  second-generation  schematic-based  design  tools  for  MEMS.  As 
shown  in  Figure  4(b),  the  atomic  elements  include  beams,  plates,  comb-drives,  and  anchors,  which 
can  be  combined  to  form  higher-level  schematics  as  in  Figure  4(a).  The  ability  to  link  elements 
from  the  catalog  at  multiple  levels  of  the  design  hierarchy  through  a  common  representation  enables 
mixed-level  simulation  of  MEMS. 

Task  5:  Technology  Integration  and  Transfer  (Microcosm) 

The  tools  and  methodologies  developed  in  this  project  were  transferred  to  Microcosm  during 

annual  workshops  held  at  Microcosm.  The  first  year  workshop  ran  though  the  entire  summer  allow¬ 
ing  close  student  interaction  [3][4][5]  and  choice  of  interfaces  between  the  student-level  tasks.  Sub¬ 
sequent  workshops  ran  for  a  single  week  and  provided  students  with  an  intensive  environment  for 
demonstrating  their  latest  developments,  as  well  as  integrating  between  selected  development  activ¬ 
ities.  Microcosm  has  licensed  the  fast  sensitivity  computation  developed  at  MIT,  and  is  licensing  the 
NODAS  schematic  representation  and  behavioral  models  and  the  synthesis  modules  from  Carnegie 
Mellon.  NODAS  will  be  packaged  with  Microcosm’s  MEMSYS  system-level  design  environment, 


and  will  complement  their  existing  macromodel  generation  features. 


3.  Deliverable  Checklist 

•  Q4-97:  2D  microresonator  synthesis  module  (CMU) 

•  Q4-97:  Semiconductor  process  representation  with  statistical  extensions  &  MUMPS  implemen¬ 
tation 

•  Q 1  -98:  3D  microresonator  synthesis  module  (CMU) 

•  Q 1  -98:  First-generation  MEMS  MAST  library  -  NODAS  1 .0  (CMU) 

•  Q 1  -98:  DC/steady-state  simulation  of  2D  beam/gap  systems  in  MATLAB  -  Sugar  0.4  (UCB) 

•  Q 1  -98:  MEMCAD  4.0  (improved  macromodeling  technology  for  use  in  MEMSYN  efforts) 
(MIT/Microcosm) 

•  Q2-98:  2D  shape  synthesis  for  specified  force-deflection  behavior  (U.  Penn) 

•  Q3-98:  Transient  analysis  planar  beam/gap  systems  in  MATLAB  -  Sugar  (UCB) 

•  Q3-98:  Nonlinear  beam  modeling  in  MATLAB  -  Sugar  (UCB) 

•  Ql-99:  Fast  geometric  sensitivity  computation  (MIT) 

•  Ql-99:  MEMCAD  4.5  (insertion  of  MEMS  schematic  technology)  (Microcosm/CMU) 

•  Ql-99:  Simulation  of  nonlinear  3D  beam/gap  systems  -  Sugar  (UCB) 

•  Q2-99:  Microaccelerometer  synthesis  module  (CMU) 

•  Q2-99:  Second-generation  analog-HDL  schematic  library  (CMU) 

•  Q3-99:  Layout-based  schematic  library  (CMU) 

•  Q3-99:  Shape  optimization  of  a  non-linear  mechanical  structure  (U.  Penn) 

•  Q4-99:  CMOS-MEMS  microaccelerometer  synthesis  module  (CMU) 
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7.  Technology  Transfer  Plans 

Microcosm  Technologies  has  licensed  some  of  the  tools  and  methodologies  developed  in 
this  project.  MEMScaP  has  expressed  interest  in  licensing  some  of  these  tools  and  methodologies. 
In  particular,  the  NODAS  schematic  representation,  as  well  as  the  schematic-centric  design  environ¬ 
ment  it  has  spawned  will  be  included  in  MEMSCAP’s  Kanaga  design  kit  for  MUMPS  and  in  a  yet  to 
be  developed  CMOS-MEMS  design  kit  for  the  DARPA  MTO  MEMS  project  on  “Application-Spe¬ 
cific  Integrated-MEMS  Process  Service.” 

8.  Future  Extensions 

To  address  the  link  between  the  schematic  and  layout  representations,  a  layout  to  schematic 
extractor  is  being  developed  under  the  joint  DARPA  MTO  MEMS  and  Composite  CAD  project  on 
“Integrated  MEMS  Inertial  Measurement  Unit,”  and  schematic-driven  layout  generation  is  being 
developed  under  the  DARPA  MTO  MEMS  project  on  “Application-Specific  Integrated  MEMS  Pro¬ 
cess  Service.” 

Extensions  of  the  hierarchical  representation  methodology  to  enable  layout  generation  and 
design  rule  checking  is  being  developed  under  the  auspices  of  the  DARPA  MTO  “Application-Spe¬ 
cific  Integrated  MEMS  Process  Service”  project,  and  MEMS  and  parasitic  extraction  is  being  devel¬ 
oped  in  the  DARPA  MTO  “Integrated  MEMS  Inertial  Measurement  Unit”  project.  These  efforts 
link  the  schematic  representation  developed  in  this  project  with  the  existing  layout  and  solid  model 
representations  allowing  the  MEMS  designer  to  choose  their  preferred  mode  of  design  entry  and 
design  simulation  independently. 
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ABSTRACT 

In  this  survey  paper  we  describe  and  contrast  three  different  approaches  for  extending  circuit  simulation  to 
include  micromachined  devices.  The  most  commonly  used  method,  that  of  using  physical  insight  to  develop 
parameterized  macromodels,  is  presented  first.  The  issues  associated  with  fitting  the  parameters  to  simulation 
data  while  incorporating  design  attribute  dependencies  is  considered.  The  numerical  model  order  reduction 
approach  to  macromodeling  is  presented  second,  and  some  of  the  issues  associated  with  fast  solvers  and  model 
reduction  are  summarized.  Lastly,  we  describe  the  recently  developed  circuit-based  approach  for  simulating 
micromachined  devices,  and  describe  the  design  hierarchy  and  the  use  of  a  catalog  of  parts. 


1  Introduction 

Decades  of  enormous  research  and  capital  investment  in  VLSI  technology  has  made  it  possible  to  put  more 
than  a  million  transistors  on  a  square  centimeter  of  silicon,  and  that  investment  is  now  also  making  it  possible 
to  fabricate  devices  with  micron-scale  moving  parts.  The  specific  techniques  used  to  fabricate  such  vanishingly 
small  moving  parts  is  often  referred  to  as  micromachining,  and  the  potential  impact  of  micromachining  is  hard 
to  understate.  Micromachined  devices  will  play  a  key  role  in  making  the  now  pervasive  computer  technology 
interact  more  directly  with  the  physical  world.  Micromachined  devices  are  already  providing  such  physical- 
computer  interfaces:  micromachined  accelerometers  are  used  in  automobile  automatic  airbag  deployment 
systems  [1],  micromachined  million  mirror  arrays  are  used  in  computer  projection  displays  [2],  and  centimeter¬ 
sized  pressure  sensors  are  used  in  a  range  of  industrial  control  applications  [3]. 

Researchers  in  almost  every  engineering  and  scientific  discipline  are  examining  ways  to  harness  the  ability 
to  fabricate,  at  low  cost,  centimeter-sized  systems  with  of  hundreds  of  thousands  of  mechanical  parts  and 
transistors.  Microresonators,  which  can  replace  bulky  passive  components  in  communication  circuits,  may 
usher  in  wristwatch-sized  cell  phones  (for  better  or  worse)  [4];  active  research  on  microfluidic  valves,  pumps 
and  mixers  may  lead  to  single-chip  chemical  analysis  systems  which  could  be  used  to  make  “in-vitro”  medical 
diagnostic  equipment  or  pocket-sized  chemical  agent  detectors  [5];  research  on  microfabricated  turbines  and 
generators  [6]  may  lead  to  an  alternative  to  batteries  for  portable  energy;  and  microfabricated  parts  small 
enough  to  capture  and  hold  individual  biological  cells  will  accelerate  progress  in  both  medical  and  scientific 
research  [7], 

Over  the  last  decade  there  has  been  extensive,  and  successful,  research  focussed  on  developing  and  exploiting 
micromachining,  though  there  are  very  few  high-volume  micromachined  products.  In  addition,  almost  all  the 
research  in  applying  micromachining  technology  has  been  carried  out  by  specialists  with  many  years  of  focussed 
training.  In  contrast,  integrated  circuit  designers  do  not  need  such  a  high  level  of  specialization.  Instead,  they 
rely  on  a  coordinated  suite  of  synthesis  and  verification  tools  that  makes  it  possible  to  design  an  application- 
specific  circuit  with  high  confidence  of  first-pass  success,  even  without  becoming  an  expert  in  semiconductor 
fabrication.  The  current  situation  for  micromachined  device  designers  is  very  different.  These  designers  must 
know  the  fabrication  process  intimately,  and  may  even  have  to  design  their  own  process.  In  addition,  the 
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Figure  1:  System-level  behavioral  model  of  a  multi-resonator  filter 


design  tools  available  are  often  limited  and  provide  only  domain-specific  simulation  or  rudimentary  layout 
editing.  The  combination  of  inadequate  computer-aided  design  tools  and  rapidly  evolving  process  technology 
has  created  an  expertise  barrier  that  excludes  nonspecialists  who  would  bring  important  application  expertise. 
Unless  this  expertise  barrier  is  lowered,  primarily  through  vastly  improved  computer-aided  design  tools,  it 
seems  unlikely  that  the  potential  of  micromachining  to  impact  so  many  different  disciplines  will  be  achieved. 

Although  the  need  for  design  tools  for  micromachining  has  been  recognized  for  well  over  a  decade,  progress 
has  been  stymied  by  a  problem  whose  difficulty  has  been  persistently  underestimated.  To  introduce  this 
problem,  consider  that  for  nearly  thirty  years,  integrated  circuit  designers  have  relied  on  circuit  simulation. 
This  one  tool  has  nearly  eliminated  the  need  to  build  prototype  circuits  in  order  to  find  major  design  flaws. 
One  reason  for  the  success  of  circuit  simulation  is  that  its  input  is  the  same  schematic  diagram  that  designers 
use  to  reason  about  the  circuit,  and  the  simulator’s  output  is  roughly  the  same  as  would  be  produced  by 
prototyping  the  circuit  and  then  measuring  all  the  voltages  and  currents.  The  problem  for  micromachined 
designers  is  that  there  is  no  equivalent  of  a  circuit  simulator,  and  no  equivalent  of  a  schematic  language  to 
describe  the  device  to  a  simulator,  if  such  a  simulator  existed.  Simulator  extension  languages  like  VHDL- 
AMS  [8]  can  greatly  simplify  the  mechanics  of  incorporating  models  for  micromachined  devices  into  circuit 
simulators,  but  they  do  not  address  a  more  fundamental  problem.  In  a  traditional  circuit  schematic,  elements 
interact  only  at  nodes,  and  the  physical  position  of  elements  has  limited  impact  on  performance.  Neither  of 
these  these  circuit-oriented  concepts  translate  directly  to  micromachined  device  design. 

The  problem  of  how  best  to  extend  circuit  simulation  to  include  micromachined  devices  is  fundamental, 
and  as  yet,  unsolved.  For  this  reason,  in  this  paper  we  will  focus  on  the  emerging  approaches  to  simulation. 
In  order  to  make  some  of  the  issues  clearer,  we  will  start  in  the  next  section  with  a  brief  description  of  a  filter 
example  which  uses  a  micromachined  device.  Then  in  section  3,  we  will  describe  the  currently  most  widely 
used  approach  to  extending  circuit  simulation,  that  of  generating  semi-analytical  macromodels  for  each  type 
of  micromachined  device.  Then,  in  section  4,  we  will  discuss  the  desirability  and  difficulties  of  replacing  the 
semi-empirical  macromodeling  approach  with  a  purely  numerical  approach  based  on  computer  simulation  and 
model-order  reduction.  In  section  5,  we  will  approach  the  simulation  problem  from  the  specification  side,  and 
discuss  a  hierarchy  of  elements  and  a  schematic  description  for  a  certain  classes  of  micromachined  devices. 
Finally,  in  our  conclusions,  we  try  to  tie  together  these  separate  approaches  and  loosely  conjecture  about 
where  the  field  is  going. 

2  Example  and  Background 

In  this  section  we  describe  a  design  example  in  order  to  help  illustrate  the  difficulties  in  developing  extensions 
to  a  circuit  simulator  for  micromachined  devices.  The  example  is  a  bandpass  filter  which  uses  a  series  of  comb- 
drive  micromachined  resonators  [9],  shown  in  a  high-level  form  in  figure  1.  The  high-level  diagram  is  best 
described  by  tracing  from  input  to  output.  The  input  n,n  in  figure  1  is  connected  to  a  triangle  which  represents 
a  transistor  amplifier.  The  parallel  plates  adjacent  to  Fc  indicates  an  electrical  to  mechanical  conversion.  The 
force  Fc  accelerates  the  first  mass  in  a  spring-coupled  cascade  of  spring-mass-dashpot  resonators.  Finally,  the 
parallel  plates  adjacent  to  C  indicates  a  mechanical  to  electrical  conversion  which  feeds  an  transistor  amplifier 
which  generates  vout. 
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Figure  2:  Overhead  view  of  the  lateral  microresonator  (Figure  thanks  to  C.  Nguyen  and  R.  Howe,  permission 
not  yet  confirmed). 


Figure  3:  SEM  of  an  integrated  CMOS  resonator  (Figure  thanks  to  C.  Nguyen  and  R.  Howe,  permission  not 
yet  confirmed). 

In  order  to  better  understand  the  filter  example,  consider  a  single  comb-drive  lateral  microresonator,  a 
layout  is  shown  in  figure  2.  An  SEM  of  the  fabricated  device  is  shown  in  figure  3.  The  polysilicon  resonator 
structure,  which  appears  white  in  the  SEM  picture  and  gray  in  the  top  view  diagram,  has  been  released  from 
the  substrate  underneath  except  at  certain  attachment  points.  The  thick  black  lines  in  figure  2  are  used  to 
show  where  the  polysilicon  structures  are  attached  to  the  substrate,  or  anchored.  As  is  clear  from  figure  2, 
the  structure  has  thred  separately  anchored  parts:  a  left  comb,  a  right  comb,  and  a  dual-comb  central  shuttle 
which  is  anchored  to  a  only  through  thin  polysilicon  beams.  The  thin  beams  serve  two  purposes.  They  act 
as  springs  and  allow  the  central  shuttle  to  oscillate  from  left  tp  right,  and  they  provide  a  conductive  path 
between  the  centred  shuttle  and  a  fixed  conducting  plate  held  at  a  bias  potential  Vp.  The  interdigited  combs 
generate  electrostatic  forces  which  pull  the  shuttle  to  the  left  when  Vi  >  Vg  and  to  the  right  when  vq  >  t>j, 
assuming  both  w,  and  vq  are  larger  than  VP.  If  out-of-phase  sinusoidally  varying  voltages  are  applied  to  V{  and 
tty,  along  with  a  dc  offset,  then  the  amplitude  of  the  central  shuttle’s  steady-state  oscillation  will  be  strongly 
frequency  dependent. 

As  the  diagram  in  figure  2  suggests,  all  that  is  needed  to  include  the  resonator  in  a  circuit  simulator  is  to 
determine  a  relationship  between  the  currents  i,  and  i0  and  the  voltages  Vp,  v{,  and  v0.  And  at  least  formally, 
the  needed  current-voltage  relationship  <;an  be  derived  by  determining  the  mechanical  material  properties  and 
then  solving  a  coupled  system  of  time-dependent  partial  differential  equations  on  a  moving  boundary.  In 
particular,  the  shuttle  accelerates  and  the  tethers  bend  elastically  in  response  to  forces  generated  by  exterior 
electric  fields  and  viscous  drag.  The  drag  will  not  be  exactly  zero  if  the  resonator  is  packaged  in  a  vacuum, 
because  there  are  still  mechanical  energy  loss  mechanisms  which  create  an  effective  drag  force. 

Though  the  statement  in  the  preceding  paragraph  is  true,  it  is  hides  many  of  the  important  difficulties. 
Determining  a  device’s  three-dimensional  structure '  and  associated  material  properties  requires  a  detailed 
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understanding  of  the  fabrication  processes  as  well  as  a  set  of  carefully  designed  experiments  [10].  Solving 
coupled  systems  of  three-dimensional  time-dependent  partial  differential  equations  with  a  complicated  moving 
boundary  requires  sophisticated  numerical  techniques  and  a  great  deal  of  computer  time  [11,  12].  Finally, 
developing  a  current-voltage  relation  for  the  resonator  in  figure  2  may  not  even  be  an  appropriate  goal.  To 
see  why,  reconsider  the  original  filter  with  the  three  stage  resonator  shown  in  figure  1.  In  the  multi-stage 
resonator,  the  single  stage  resonators  are  coupled  together  by  springs  which  are  implemented  using  thin 
polysilicon  tethers.  For  the  multi-stage  resonator,  the  most  important  aspect  to  model  well  is  the  transfer 
function  from  the  input  of  the  first  resonator  to  the  output  of  the  last  resonator.  However,  there  will  be  no  way 
to  arrive  at  that  transfer  function  by  “composing”  the  previously  mentioned  single  resonator  models.  Instead, 
an  entirely  new  model  will  be  needed  for  a  two-stage  resonator,  and  then  another  model  for  a  three-stage 
resonator,  and  yet  another  model  for  a  four  stage  resonator,  et  cetera.  And  if  these  models  are  going  to  be 
derived  by  solving  time-dependent  partial  differential  equations  for  structures  as  geometrically  complicated  as 
a  multistage  resonator,  the  computer  time  required  may  cast  a  more  positive  light  on  building  prototypes. 

In  order  to  assess  the  importance  of  issues  like  deriving  structure  and  material  properties  from  layout  and 
process  information,  the  computational  cost  of  partial  differential  equation  solution,  or  model  composibility,  it 
is  worth  recalling  that  for  integrated  circuit  design,  simulator  use  can  be  divided  into  two  broad  classes.  Early  in 
the  circuit  design  process,  during  a  synthesis  or  optimization  phase,  many  alternatives  are  being  considered,  and 
designs  are  typically  represented  only  with  a  schematic.  That  is,  circuit  element  interconnection  is  specified, 
but  no  layout  information  exists.  The  simplicity  of  the  schematic  representation  both  builds  intuition  and 
accelerates  examining  alternatives,  though  layout  parasitics  are  either  ignored  or  crudely  estimated.  As  the 
design  matures,  when  the  circuit  layout  has  been  determined,  a  verification  phase,  begins.  Circuit  simulators 
are  then  combined  with  layout  extraction  techniques  to  check  that  the  final  layout  results  in  a.  circuit  with  the 
desired  performance.  Such  a  two  stage  approach  also  seems  to  be  a  natural  fit  to  designing  the  microresonator 
filter.  It  would  be  very  efficient  if  most  of  the  layout  details  could  be  avoided  while  examining  alternatives 
such  as:  fewer  or  more  resonator  stages;  fewer  or  more  comb  fingers;  heavier  or  lighter  shuttles;  and  longer 
or  shorter  or  more  serpentine  tethers.  Then  only  during  the  verification  phase  would  it  be  necessary  to  work 
with  the  layout  and  combine  extraction  techniques  with  simulation. 

3  Semi- Analytical  Macromodeling 

By  far  the  most  common  approach  to  including  a  micromachined  device  in  circuit  simulation  is  to  analyze  the 
device  approximately,  so  as  to  generate  a  macromodel  in  the  form  of  either  a  circuit  or  a  low-order  system  of 
differential  equations  [13].  Generating  the  form  of  these  models  requires  the  device  designer’s  physical  insight, 
and  can  be  as  much  art  as  science,  though  issues  such  as  energy  conservation  can  provide  guidelines  [14]. 
Once  the  form  of  the  macromodel  has  been  discerned,  then  values  for  the  parameters  must  be  determined. 
These  parameters  can  be  determined  analytically,  or  by  experiment,  or  by  using  numerical  simulation.  The 
decomposition  between  macromodel  form  and  parameterization  is  not  a  precise  one,  and  is  best  described  by 
example.  Below,  a  simple  macromodel  form  for  the  single  resonator  example  of  figure  2  is  derived,  and  then 
several  alternatives  for  determining  model  parameters  are  examined.  The  merits  and  deficits  of  semi-analytic 
macromodels  are  then  described  in  general. 

3.1  Example  Model  Form 

In  order  to  develop  a  model  for  the  resonator  which  can  be  incorporated  in  a  circuit  simulator,  the  resonator 
voltages  must  be  related  to  the  resonator  currents.  Resonators  are  usually  modeled  as  RLC  circuits  [9]  as  such 
models  help  develop  intuition.  A  differential  equation  model  is  developed  below  because  the  setting  is  more 
generally  applicable.  To  begin,  from  figure  2,  the  currents  i,  and  i0  can  be  related  to  the  voltages  v,,  v0  and 
Vp  by  first  noting  that 

ii(t)  =  ^Ql(v* -Vp,x)  i0(t)=  ^Qr(vo-Vp,x)  (1) 

where  x  is  the  displacement  from  center  of  the  dual-comb  shuttle,  and  Qr  and  Qi,  are  the  net  charges  on  the 
left  and  right  anchored  combs,  respectively.  A  simple  parallel  plate  analysis  suggests  that  the  comb  capacitance 
is  an  affine  function  of  the  displacement  x,  in  which  case  the  comb  charges  will  satisfy  an  equation  of  the  form 

QdVi  -  Vp,x)  «  (Co  -  C\x)(vi  -  Vp)  Qr(v,  -  Vp,x )  «  (Co  +  Cvx)(v0  -  Vp)  (2) 

where  Co  is  the  comb- pair  capacitance  at  a:  =  0  and  C\  is  the  .r-derivative  of  that  capacitance.  Note  that  there 
is  only  one  Co  and  one  Ci  so  we  have  exploited  the  left-right  symmetry  in  the  problem.  Finally,  a  very  simple 
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spring-mass  mechanical  model  for  the  comb  suggests  that  the  shuttle  displacement,  x,  satisfies  an  equation  of 
the  form 

M^x  +  Kdjtx  +  K.x  =  Ke  ((Vo  -  Vp)2  -  (Vi  -  Vpf)  (3) 

where  M  is  the  mass  of  the  shuttle,  Kd  is  the  drag  force  on  the  comb  generated  by  the  surrounding  fluid 
(typically  air  at  room  pressure),  K,  is  the  spring  constant  associated  with  the  thin  teathers,  and  Ke  is 
constant  which  relates  the  electrostatic  force  generated  by  the  comb  to  the  square  of  the  applied  voltage. 

3.2  Determining  Model  Parameters 

A  very  simple  analysis  of  the  microresonator  was  used  above  to  develop  a  differential  equation  system 
macromodel.  The  model  is  given  by  the  combination  of  (1),  (2),  and  (3).  It  is  worth  noting  that  the  model 
is  nonlinear  and  has  quite  a  few  parameters.  Until  the  parameters  are  set,  there  is  only  a  “form”  for  the 
macromodel.  In  the  above  example,  and  for  macromodels  in  general,  specifying  the  macromodel  form  usually 
implies:  assigning  a  set  of  state  variables,  determining  which  time  derivatives  appear,  representing  which  state 
variables  interact,  and  specifying  where  the  parameters  appear.  It  is  also  quite  common  to  include  certain 
expected  nonlinearities,  as  was  done  with  the  squared  potential  in  (3). 

The  above  macromodel  has  many  parameters,  Co,  C\ ,  M,  Kd,  Ks  and  Ke,  so  it  is  tempting  to  suggest  that 
the  model  could  fit  anything.  Since  the  macromodel  is  intended  only  as  an  example,  we  will  consider  the  issue 
of  how  to  determine  the  parameters  rather  than  focusing  on  how  to  improve  the  model.  There  are  two  main 
issues  associated  with  macromodel  parameter  selection: 

1.  Will  the  parameters  be  determined  by  physical  analysis  or  through  fitting  to  measured  or  simulated  data? 

2.  Will  a  new  set  of  parameters  be  determined  every  time  a  change  is  made  in  the  device  geometry ,  or  will 
the  parameters  be  given  as  an  explicit  function  of  design  attributes? 

Most  macromodels  use  some  combination  of  analysis  and  data  fitting  to  determine  the  parameters.  For 
example,  the  shuttle  mass  of  the  resonator,  M,  is  easily  determined  from  the  geometry  and  material  properties, 
but  would  be  difficult  to  measure  directly.  There  are  techniques  for  estimating  shuttle  drag  [15],  Kd,  though 
recent  studies  suggest  that  numerically  solving  the  Stoke ’s  equation  yields  higher  accuracy  [16].  Finite-element 
analysis  or  measurements  might  also  be  superior  to  trying  to  use  linear  beam  theory  when  attempting  to 
determine  the  spring  coefficient  KB  [10,  17].  In  general,  as  software  for  solving  partial  differential  equations 
improves,  parameter  estimation  will  be  more  heavily  based  on  results  from  simulation  rather  analytical 
techniques. 

There  are  many  aspects  of  a  microresonator  that  a  designer  can  alter  to  try  to  improve  performance 
including:  the  number  and  length  of  comb  fingers,  the  tether  lengths  and  widths,  and  the  shuttle  proof  mass. 
One  advantage  of  using  physical  analysis  to  determine  macromodel  parameters  is  that  the  analysis  usually 
reveals  an  explicit  form  for  the  dependence  of  the  parameters  on  design  attributes.  Macromodels  whose 
parameters  are  given  as  explicit  functions  of  design  attributes  are  of  obvious  value  during  the  synthesis  and 
optimization  phase  of  design.  If  the  electrostatic  force  constant,  Ke,  is  estimated  analytically  using  a  parallel- 
plate  formula,  then  the  resulting  formula  for  the  parameter  Ke  will  include  a  term  which  grows  linearly  with 
the  number  of  comb  fingers. 

Deriving  macromodel  parameters  by  fitting  to  measured  data  or  simulation  results  does  not  preclude 
generating  macromodels  whose  parameters  are  given  as  explicit  functions  of  the  design  attributes.  It  is  possible 
to  use  a  multivariate  polynomial  fitting  procedure  to  generate  these  explicit  functions,  but  the  procedure  is  not 
completely  automatic  and  requires  expert  input  [17].  To  understand  the  difficulty,  consider  a  micromachined 
device  whose  macromodel  has  a  parameter  P  that  is  dependent  on  the  value  of  d  design  attributes.  We 
will  denote  the  values  of  the  d  design  attributes  as  a  d-length  vector  u.  Then,  our  problem  becomes  one  of 
determining  an  explicit  representation  of  P(u). 

A  seemingly  straight-forward  approach  to  finding  an  explicit  representation  of  P(u)  is  to  use  a  multivariate 
polynomial.  To  see  the  difficulty  generated  by  such  an  approach,  consider  all  the  terms  of  a  second-order 
polynomial  in  three  design  attributes, 

P(m)  =  Po  +  PiUi  +  P2U2  +  P3U3  +  pdu\Ui  +  /85M1M2  (4) 

+PqU]U3  +  PjU2U2  +  PsUzUz  +  P9U3U3, 

where  u  =  [?ti  U2U3] T  are  the  design  attributes  and  the  pi’s  are  the  unknown  coefficients  of  the  multivariate 
polynomial.  As  should  be  clear  from  the  above  example,  the  number  of  terms  in  a  qth  order  d-dimensional 
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polynomial  is  proportional  to  dq .  This  implies  that  it  will  be  computationally  hopeless  to  use  multivariate 
polynomials  directly  to  represent  parameter  variation  when  the  number  of  design  attributes  exceeds  a  half 
dozen.  Instead,  the  polynomials  will  have  to  be  modified  by  “pruning”  unnecessary  terms.  Determining  which 
terms  can  be  safely  discarded  requires  significant  mathematical  and  physical  insight. 

There  is  a  second  issue  associated  with  fitting  P(u)  with  a  polynomial,  and  this  issue  is  sometimes  referred 
to  as  the  “design-of-experiments”  problem.  Consider  again  the  problem  of  fitting  P(u)  with  a  second-order 
polynomial  in  three  design  attributes.  In  order  to  compute  the  ten  unknown  polynomial  coefficients  /5« ,  -  -  ■ ,  , 

values  for  P(u)  must  be  computed  for  at  least  10  different  values  of  the  3-length  vector  u.  There  are  several 
approaches  to  determining  the  “best”  test  values  for  u  [17]  based  on  statistical  arguments,  but  the  key  difficulty 
and  its  cure  can  be  seen  by  examining  the  matrix  equation  associated  with  the  fitting.  Again  for  our  second- 
order  example,  assume  there  are  l  +  1  measurements,  l  >  9,  and  let  the  measurements  (either  real  or  from 
simulation)  be  at  points  u'  for  i  =  0  to  l.  The  matrix  equation  for  the  p's  is  then  given  by 


P(u°)  ' 
P(u') 

1  u®  u®  . 
1  u\  u\ 

■  »3U3  " 

U3U3 

'  Po  ‘ 
Pi 

P(u‘)  . 

1  ?l[  Ul2  ■ 

■  U3U3  . 

.  ft-  . 

(5) 


Note  that  the  system  will  be  square  when  the  number  of  measurements  equals  the  number  coefficients,  though 
typically  the  number  of  measurements  far  exceeds  the  number  of  coefficients. 

As  is  clear  from  examining  (5),  the  points  uq  through  «/  should  be  chosen  to  make  the  rows  (or  the  columns) 
of  the  matrix  in  (5)  as  close  to  mutually  orthogonal  as  possible.  Finding  values  of  m”s  which  generate  a  nearly 
orthogonal  matrix  can  be  accomplished  using  a  one  test  point  at  a  time  algorithm. 


3.3  Merits  and  deficits 

The  semi-analytic  approach  to  macromodeling  is  in  far  wider  use  than  the  methods  to  be  described  below. 
And,  since  this  method  is  “free-form”,  there  are  no  restrictions  as  to  what  kind  of  micromachined  devices  can 
be  modeled.  In  addition,  if  such  macromodels  are  carefully  parameterized,  they  can  be  used  to  excellent  effect 
during  the  synthesis  and  optimization  phase  of  design. 

There  are  two  difficulties  with  the  semi-analytic  macromodeling  approach.  The  most  obvious  problem  is 
that  there  is  no  standard  method  for  generating  these  macromodels,  and  the  only  way  to  determine  when 
the  models  are  sufficiently  accurate  is  by  comparing  the  macromodel’s  results  to  those  of  experiments  or  very 
detailed  simulation.  The  second  problem  is  simply  that  the  macromodeling  approach  provides  very  little 
verification.  One  can  not  “extract”  the  macromodel  from  layout,  or  add  in  parasitics.  In  addition,  since  the 
device  designer  usually  generates  the  macromodel,  there  is  no  independent  check  on  whether  an  important 
interaction  is  being  ignored. 


4  Numerical  Macromodeling 

When  the  first-  and  second-order  behaviors  of  a  micromachined  device  are  well-understood,  the  most  efficient 
strategy  for  including  the  device  in  a  circuit  simulator  is  to  develop  the  kind  of  semi-analytic  macromodel 
discussed  in  the  previous  section.  Given  how  rapidly  the  micromachining  technology  is  changing,  it  is 
rarely  possible  to  wait  for  such  device  expertise  to  develop.  And  since  it  is  almost  impossible  to  design 
systems  which  use  micromachined  devices  without  access  to  reliably  accurate  macromodels,  slow  macromodel 
development  translates  into  slow  technology  deployment.  For  this  reason,  there  has  been  a  steady  effort  over 
the  last  decade  to  develop  nearly  automatic  approaches  for  generating  accurate  macromodels  of  micromachined 
devices  starting  from  only  layout  and  process  descriptions.  Most  efforts  is  this  area  are  following  a  three  step 
approach:  [18,  11,  19] 

•  Use  modified  extrusion  to  generate  an  approximate  3-D  structure  from  a  layout  and  process  descrip¬ 
tion  [10,  20]. 

•  Use  fast  coupled-domain  3-D  simulation  techniques  to  analyze  the  entire  micromachined  device  [12,  21]. 

•  Use  a  projection-based  model-order  reduction  strategy  to  generate  macromodels  from  3-D  simulation  [22]. 

Below  we  describe  some  of  the  recent  developments  and  persistent  challenges  in  fast  coupled-domain 
simulation  and  model-order  reduction. 
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4.1  Fast  Coupled-domain  3-D  Solvers 

The  above  strategy  for  generating  macromodels  by  performing  projection-based  model-order  reduction  relies 
on  the  ability  to  simulate  an  entire  micromachined  device  in  a  reasonable  period  of  time.  In  order  to  simulate 
the  microresonator  in  figure  2,  for  example,  it  is  necessary  to  solve  a  complicated  three-dimensional  moving 
boundary  problem  which  couples  elastic,  fluidic  and  electrostatic  forces.  Simulating  such  problems  with 
standard  finite-element  methods  is  nearly  intractable,  because  it  is  necessary  to  discretize  the  volume  in  both 
the  interior  and  exterior  of  the  resonator.  In  order  to  simulate  entire  micromachined  devices,  it  was  first 
necessary  to  develop  faster  techniques  for  analyzing  the  exterior  field  problems.  Then,  it  became  necessary  to 
develop  robust  approaches  for  coupling  those  fast  techiques  to  standard  finite-element  algorithms  for  computing 
elastic  deformation. 


4.1.1  Fast  Field  Solvers 


For  most  micromachined  devices,  the  electrostatic  and  fluidic  forces  in  the  exterior  of  the  device  satisfy 
linear  partial  differential  equations.  Specifically,  the  surface  electrostatic  forces  can  be  determined  by  solving 
an  exterior  Laplace’s  equation  and  the  fluid  drag  forces  can  be  determined  by  solving  an  exterior  Stoke’s 
equation.  For  both  equations,  it  is  possible  to  derive  integral  formulations  which  avoid  the  exterior  volume 
entirely  and  instead  relates  potentials  to  forces  in  the  electrostatic  case,  and  velocities  to  forces  in  the  fluid 
case. 

More  specifically,  the  electrostatic  potential  and  the  fluid  velocity,  assuming  Stoke’s  flow,  both  satisfy  an 
integral  equation  over  the  device  surface  given  by  Green’s  theorem: 


«(; r) 
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surfaces 


G(x,  x') 


du(x') 

dn 
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9G(x,x') 
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where  u  is  either  the  electrostatic  potential  or  the  fluid  velocity,  x  is  a  point  on  the  surface,  and  j~n  is  the 
derivative  in  the  direction  normal  to  the  polysilicon  surface. 

Discretization  of  the  above  integral  equation  leads  to  a  dense  system  of  equations  which  becomes 
prohibitively  expensive  to  form  and  solve  for  complicated  problems.  To  see  this,  consider  the  electrostatics 
problem  of  determining  the  surface  charge  given  the  potential  on  conductors.  A  simple  discretization  for  the 
electrostatics  problem  is  to  divide  the  polysilicon  surfaces  into  n  flat  panels  over  which  the  charge  density  is 
assumed  constant.  A  system  of  equations  for  the  panel  charges  is  then  derived  by  insisting  that  the  correct 
potential  be  generated  at  a  set  of  n  test,  or  collocation,  points.  The  discretized  system  is  then 


Pq  =  $  (7) 

where  q  is  the  7i-length  vector  of  panel  charges,  $  is  the  n-length  vector  of  known  centroid  potentials.  Since 
the  Green’s  function  for  electrostatics  is  the  reciprocal  of  the  separation  distance  between  x  and  x' , 
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and  therefore  every  entry  in  P  is  nonzero. 

If  direct  factorization  is  used  to  solve  (7),  then  the  memory  required  to  store  the  matrix  will  grow  like 
n2  and  the  matrix  solve  time  will  increase  like  n3.  If  instead,  a  preconditioned  Krylov-subspace  method  like 
GMRES  [23]  is  used  to  solve  (7),  then  it  is  possible  to  reduce  the  solve  time  to  order  n2  but  the  memory 
requirement  will  not  decrease. 

In  order  to  develop  algorithms  that  use  memory  and  time  that  grows  more  slowly  with  problem  size,  it  is 
essential  not  to  form  the  matrix  explicitly.  Instead,  one  can  exploit  the  fact  that  Krylov-subspace  methods 
for  solving  systems  of  equations  only  require  matrix-vector  products  and  not  an  explicit  representation  of  the 
matrix.  For  example,  note  that  for  P  in  (7),  computing  Pq  is  equivalent  to  computing  n  potentials  due  to 
n  charged  panels  and  this  can  be  accomplished  approximately  in  nearly  order  n  operations  [24,  25].  To  see 
how  to  perform  such  a  reduction  in  cost,  consider  Figure  (4).  The  short-range  interaction  between  close-by 
panels  must  be  computed  directly,  but  the  interaction  between  the  cluster  of  panels  and  distant  panels  can  be 
approximated.  In  particular,  as  Figure  (4)  shows,  the  distant  interaction  can  be  computed  by  summing  the 
clustered  panel  charges  into  a  single  multipole  expansion  (denoted  by  M  in  the  figure),  and  then  the  multipole 
expansion  can  be  used  to  evaluate  distant  potentials. 

Several  researchers  simultaneously  observed  the  powerful  combination  of  integral  equation  approaches, 
Krylov-subspace  matrix  solution  algorithms,  and  fast  matrix- vector  products  [26,  32].  Perhaps  the  first 
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Figure  4:  A  Cluster  of  collocation  points  separated  from  a  cluster  of  panels 


practical  use  of  such  methods  combined  the  fast  multipole  algorithms  for  charged  particle  computations  with 
the  above  simple  discretization  scheme  to  compute  3-D  capacitance  and  electrostatic  forces  [27].  Higher-order 
elements  and  improved  efficiency  for  higher  accuracy  have  been  the  recent  developments  [21,  31].  The  many 
different  physical  domains  involved  in  micromachined  devices  has  focussed  attention  on  fast  techniques  which 
are  Green’s  function  independent,  such  as  the  precorrect-FFT  schemes  [25,  29]. 

4.1.2  Coupled- Domain  Simulation 

Self-consistent  electromechanical  analysis  of  micromachined  polysilicon  devices  typically  involves  determining 
mechanical  displacements  which  balance  elastic  forces  in  the  polysilicon  with  electrostatic  pressure  forces  on 
polysilicon  surface.  The  technique  of  choice  for  determining  elastic  forces  in  the  polysilicon  is  to  use  finite- 
element  methods  to  generate  a  nonlinear  system  equations  of  the  form 

F(u)-P(t»,fl)  =  0.  (9) 

where  u  is  a  vector  of  finite-element  node  displacements,  F  relates  node  displacements  to  stresses,  and  P  is 
the  force  produced  by  the  vector  representing  the  discretized  surface  charge  q.  Note  that  as  the  structure 
deforms,  the  pressure  changes  direction,  so  P  is  also  a  function  of  u.  One  can  view  this  mechanical  analysis 
as  a  “black  box”  which  takes  an  input,  q,  and  produces  an  output  u  as  in 

«  =  HM(q).  (10) 

In  order  to  determine  the  charge  density  on  the  polysilicon  surface  due  to  a  set  of  applied  voltages,  one  can 
use  a  fast  solver,  as  described  above.  One  can  view  the  electrostatic  analysis  as  a  “black  box”  which  takes,  as 
input,  geometric  displacements,  it,  and  produces,  as  output,  a  vector  of  discretized  surface  charges,  q,  as  in 

q  =  He(u).  (11) 

Self-consistent  analysis  is  then  to  find  a  it  and  q  which  satisfies  both  (10)  and  (11). 

A  simple  relaxation  approach  to  determining  a  self-consistent  solution  to  (10)  and  (11)  is  to  successively 
use  (10)  to  update  displacements  and  then  using  (11)  to  update  charge.  Applying  (10)  implies  solving  the 
nonlinear  equation,  (9),  typically  using  Newton’s  method. 

Although  the  relaxation  method  is  simple  it  often  does  not  converge.  Instead  one  can  apply  Newton’s 
method  to  the  system  of  equations 
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in  which  case  the  updates  to  charge  and  displacement  are  given  by  solving 
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The  above  method  is  referred  to  as  a  multi-level  Newton  method  because  forming  the  right-hand  side  in  (13) 
involves  using  Newton’s  method  to  apply  Hm- 

In  order  to  solve  (13),  one  can  apply  a  Krylov-subspace  iterative  method  such  as  GMRES.  The  important 
aspect  of  GMRES  is  that  an  explicit  representation  of  the  matrix  is  not  required,  only  the  ability  to  perform 
matrix-vector  products.  As  is  clear  from  examining  (13),  to  compute  these  products  one  need  only  compute 
A q  and  ^^-Au.  These  products  can  be  approximated  by  finite  differences  as  in 

9Hm  Aq  Hm('Q  +  aAq^  ~  Hm^  Q41 

dq  ~  a  '  ' 

where  a  is  a  very  small  number.  Therefore,  this  matrix-free  multilevel-Newton  method  [34]  can  treat  the 
individual  solvers  as  black  boxes.  The  black  box  solvers  are  called  once  in  the  outer  Newton  loop  to  compute 
the  right  hand  side  in  (13)  and  then  called  once  per  each  GMRES  iteration.  Computing  Hm{q  +  &d\ )  means 
using  an  inner  loop  Newton  method  to  solve  (9),  which  is  expensive,  though  improvements  can  be  made  [12]. 
An  important  advantage  of  matrix-free  multilevel-Newton  methods  is  that  it  is  not  necessary  to  modify  either 
the  mechanical  or  electrostatic  analysis  programs. 

4.1.3  Analyzing  a  Micromachined  Resonator 

In  order  to  demonstrate  that  the  above  techniques  make  it  computationally  feasible  to  simulate  an  entire 
resonator,  we  now  present  results  from  a  multilevel-Newton  coupled  electromechanical  solver  [30].  The  program 
uses  the  precorrected-FFT  accelerated  integral  equation  solver  [29]  with  planar  triangular  panels  to  compute 
the  electrostatic  forces.  A  finite-element,  mixed  rigid/elastic  mechanical  analysis  program  using  20  noded 
isoparametric  brick  elements  [33]  is  used  to  compute  displacments.  The  multilevel-Newton  method  uses  to 
pressure  sensitivites  to  improve  efficiencies. 

An  eighteen  finger  polysilicon  resonator  is  shown  in  Figure  5.  In  this  resonator,  the  central  shuttle  is  is 
suspended  by  400  micron  long  folded  beams  with  a  uniform  thickness  of  1.94  microns  and  finger  dimensions 
of  13.8x4.6  microns.  The  central  shuttle  and  an  underside  fixed  plate  (not  shown)  were  set  to  zero  volts,  and 
a  drive  voltage  was  applied  to  the  right-  and  left-hand  side  combs  (also  not  shown).  For  this  example,  the 
Young’s  modulus  of  the  polysilicon  was  determined  to  be  150 MPa,  and  the  poisson  ratio  was  0.3. 

The  effect  of  varying  the  separation  of  the  suspension  beams,  shown  as  L  in  Figure  5,  on  levitation  was 
investigated  using  the  coupled-domain  solver.  The  results  are  plotted  in  Figure  6,  which  shows  levitation 
(motion  normal  to  the  substrate)  as  a  function  of  applied  comb  drive  voltage.  Levitation  in  resonators  is 
to  be  avoided,  because  raising  the  central  shuttle  causes  a  misalignment  of  the  interdigitated  fingers.  This 
misalignment  reduces  the  resulting  electrostatic  forces,  and  may  also  allow  the  central  shuttle  to  twist  and 
collided  with  the  side  combs.  The  simulation  results  plotted  in  Figure  6  show  that  levitation  in  the  resonator 
is  can  be  nearly  as  large  as  the  resonator  thickness,  and  that  changing  separation  L  of  the  suspension  beam 
inperceptibly  effects  levitation.  The  simulation  was  run  on  a  Sun  Ultra  30,  and  each  load  step  required  70 
minutes  of  CPU  time. 


4.2  Model-Order  Reduction 

Many  micromachined  devices  are  nonlinear,  and  extracting  dynamically  accurate  nonlinear  macromodels  from 
simulation  is  a  relatively  open  problem.  For  this  reason,  there  has  been  much  current  interest  in  developing 
nonlinear  model-order  reduction  strategies  [35,  36,  37],  To  better  describe  the  challenges  in  nonlinear  model- 
reduction,  consider  simulating  the  dynamics  of  a  fixed-fixed  beam  in  a  fluid  (air).  Figure  7  [38]  shows  the 
front  view  of  the  structure.  When  a  voltage  is  applied,  the  flexible  top  plate  deforms  downward  due  to  the 
generated  electrostatic  force,  and  the  squeezed  air  in  the  gap  damps  the  plate  motion  through  a  back  pressure 
force.  The  exact  deformation  of  the  top  plate  due  to  the  applied  voltage  is  sensitive  to  the  ambient  pressure 
of  the  air,  so  this  structure  can  be  used  as  a  pressure  sensor. 

Following  Hung  et  al  [38],  the  dynamic  behavior  of  this  coupled  electro-mechanical-fluid  system  can  be 
modeled  with  the  ID  Euler  beam  equation  (15)  and  the  2D  Reynold’s  squeeze  film  damping  equation  (16): 
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Figure  6:  Levitation 


2  um  of  poly  Si  r(t)  middle  point 

0.5  um  poly  Si  displacement  as  output 


0.5  um  SiN 


2.3  um  gap 
filled  with  air 


Figure  7:  The  fixed-fixed  beam  structure.  The  x  and  z  axes  are  parallel  to  the  length  and  width  of  the  beam, 
respectively,  and  the  y  axis  points  into  the  page. 


where  x,  y  and  z  are  as  shown  in  Figure  7,  E  is  the  Young’s  modulus,  I  is  the  moment  of  inertia,  S  is  the  stress 
coefficient,  p  is  the  density,  pa  is  the  ambient  pressure,  p  is  the  viscosity,  K  is  the  Knudsen  number,  w  is  the 
width  of  the  beam  in  the  y  direction,  u  —  u(x,  t)  is  the  height  of  the  beam  above  the  substrate,  and  p(x,  y,  t)  is 
the  pressure  distribution  in  the  fluid.  Finally,  the  electrostatic  force  is  approximated  assuming  nearly  parallel 
plates  and  is  given  by  Fetec  “  where  V  is  the  applied  voltage. 

Spatial  discretization  of  (15)  and  (16)  leads  to  a  large  nonlinear  system  of  the  form 

x  =  f(x)  +  bTu(t)  y  =  cTx  (17) 

were  x  is  an  n-length  state  vector,  in  this  case  the  vector  of  displacements  and  their  time-derivatives.  The 
function  f,  which  maps  an  n-length  vector  to  an  n-length  vector,  represents  the  spatially  discretized  partial 
differential  equation.  The  above  system  with  a  nonlinear  state  equation  is  referred  to  as  the  “original”  system 
which  will  be  reduced  to  a  much  smaller  system.  The  applied  voltage  generates  u(t),  the  input  of  the  system. 
The  output  of  the  system  is  y(t),  and  is  chosen  to  be  the  beam’s  center  point  displacement. 

4.2.1  Numerical  Model  Reduction 

The  goal  of  numerical  model-order  reduction  is  to  generate  a  model  with  many  fewer  than  n  states  which  still 
preserves  the  input/output  behavior  of  the  original  system.  Almost  all  the  numerical  model-order  reduction 
strategies  are  based  on  a  change  of  variables, 

F^x^  —  x,  (18) 

where  xq  is  a  q-length  vector  (q  is  assumed  much  much  less  than  n),  and  Vq  is  n  x  q  orthonormal  matrix  whose 
columns  represent  important  “modes”.  Then,  the  matrix  Vq  represents  a  transformation  from  the  original 
to  the  reduced  coordinate  system.  Substituting  the  change  of  variables  in  (17)  and  multiplying  the  resulting 
equation  by  Vj  yields 

**  =  VJ fiYqXq)  +  Vqlm{t)  y  =  cTVqxq.  (19) 

It  should  be  noted  that  the  dynamical  system  could  have  been  multiplied  by  a  second  transformation  matrix, 
Uq,  leading  to  a  wider  range  of  algorithms. 

Buried  in  (19)  are  the  two  key  model-order  reduction  issues.  First,  one  must  select  a  good  change  of 
variables  so  that  the  input/output  behavior  is  captured  by  the  q  states  in  the  reduced  system.  Second,  and 
perhaps  less  obvious,  one  must  have  a  representation  of  Vj f{Vq-)  that  can  be  efficiently  stored  and  evaluated. 
For  example,  suppose  n  =  100,000  and  q  =  10.  Then  computing  Vj f(Vqxq)  explicitly  would  require  on  the 
order  of  100, 000  operations,  and  that  hardly  satisfies  the  efficiency  goal  of  model-order  reduction.  If  /  were 
linear,  so  that  /(x)  =  Ax  where  A  is  an  n  x  n  matrix,  then  the  representation  problem  is  easily  solved.  To 
see  this,  consider  that  for  the  linear  case 

fq(xq)  =  Vqrf(Vqxq)  =  Vq  AVqxq  =  Aqxq,  (20) 

where  fq(-)  is  a  function  which  maps  a  ^-length  vector  to  a  q-length  vector  and  denotes  the  general  nonlinear 
representation  of  the  reduced  model  in  the  reduced  variables.  The  matrix  Aq  is  a  representation  of  the  reduced 
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model  which  can  be  used  when  the  problem  is  linear.  As  Is  clear  from  the  equations,  Aq  is  an  easily  computed 
q  x  q  matrix.  For  the  example  numbers  above,  using  the  Aq  matrix  representation  to  compute  V'qr f(Vqxq) 
costs  only  100  operations  instead  of  order  100,000. 

Returning  to  the  first  issue,  selecting  the  change  of  variables  or  equivalently  chosing  Vq,  there  are  a 
number  of  methods.  If  the  problem  is  linear,  the  methods  for  determining  Vq  include:  examining  Krylov- 
subspaces  [39,  22,  40],  selecting  from  orthogonalized  time-series  data  [35],  or  computing  singular  vectors  of  the 
underlying  differential  equation  Hankel  operator  [41].  The  approach  based  on  using  time-series  data  extends 
directly  to  the  nonlinear  cases,  and  the  Krylov-subspace  and  Hankel  operator  approaches  can  be  extended 
to  the  nonlinear  case  by  linearizing  the  system  only  for  the  purpose  of  computing  Vq,  and  then  applying 
the  change  of  variables  to  the  original  nonlinear  system.  As  shown  in  [37],  linearization  approaches  can  be 
ineffective  and  better  strategies  may  exist. 

Regardless  of  how  the  V9’s  are  computed,  for  nonlinear  problems  there  is  still  the  difficulty  of  finding  an 
efficient  representation  for  Vq  f(Vq).  One  approach  is  to  assume  the  reduced  model  is  a  multidimensional 
quadratic  [43,  42],  in  which  case 

VtTf(Vqxq)  =  fq(xq)  *  J^xq  +  xTqJ™xq  (21) 

where  Jq^  is  the  q  x  q  Jacobian,  or  first  derivative,  of  fq  and  Jq 2)  is  a  q  x  q  x  q  second  derivative  of  fq. 
Both  Jq}i  and  Jq2^  are  easily  computed  from  /  by  finite-differences,  though  q2  function  evaluations  are  needed 
to  evaluate  Jq 2)  [42].  If  higher  order  nonlinearities  are  required,  such  as  cubic  or  quartic  terms,  the  above 
strategy  becomes  computationally  ineffective.  The  difficulty  arises  from  the  fact  that  there  are  qk+y  entries  in 
the  kth  derivative  of  fq,  so  generating  a  tenth-order  reduced-order  model  and  including  all  the  quartic  terms 
requires  a  representation  with  over  100,000  entries.  It  is  possible  to  use  heuristics  to  prune  the  higher-order 
nonlinearities,  so  that  only  a  small  fraction  of  the  qk+'  terms  are  retained.  Equivalently,  the  problem  is  one 
of  determining  a  sparse  representation  for  Jqk\ 

An  alternative  view  of  the  nonlinear  model  reduction  problem  can  be  developed  for  the  case  where  the 
original  nonlinear  function,  /,  can  be  represented  as  the  gradient  of  a  scalar  function  [45].  That  is, 

f(x)  =  Vx4>(x)  (22) 


where  4>  maps  an  ?i-length  vector  to  a  scalar.  Such  representations  occur  naturely  for  second-order  energy- 
conserving  systems, 

£  =  Vx4>(x)  +  bTu(t),  (23) 

where  <t>(x)  is  derived  by  constructing  the  system’s  associated  Hamiltonian.  For  such  systems,  the 
representation  problem  can  be  reduced  to  a  single  fitting  problem  by  noting  that 

f„(xq)  =  V'rf(Vqxq)  =  VjV X4>(yvxq)  =  VXq4>q{xq),  (24) 

where  (j>q  maps  an  (/-length  vector  to  a  scalar.  Then,  the  scalar  function  of  q  variables,  <j>q,  can  be  approximately 
represented  using  a  (/-dimensional  kth- order  polynomial 

Hi  H* 

4>q(xq)  as  ^2  •••  52  (25) 

>1=0  i»=0 


If  one  represents  /  directly  using  derivatives,  as  in  the  previous  paragraph,  there  are  order  qk+l  terms  in  the 
reduced  model.  Using  a  polynomial  to  represent  4>  up  to  order  k  requires  only  order  qk  terms,  so  it  would 
seem  that  exploiting  /  =  V4>  results  in  a  saving  of  only  a  factor  of  q.  However,  one  can  fit  <p  with  a  rational 
function 


SqLo  E>k=o  a>i .■■■>* xq\ 11  ••• x • 


<t>q{xq) 


ti  =o  ***  2-^1* 


Hk 
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(26) 


and  rational  function  representations  can  be  effectively  much  higher  order  than  k  without  the  commersurate 
increase  in  cost  [45] . 


4.2.2  Clamped  beam  example 

We  now  present  the  results  of  comparing  the  reduced-order  models  generated  using  the  Arnoldi  method  and  a 
finite  difference  solution  of  the  original  non-linear  governing  equations  which  is  provided  by  Hung  [38].  Hung 
verified  this  non-linear  solution  with  experimented  data. 
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Figure  8:  A  comparison  of  the  time  responses  of  the  non-linear  model,  linear  model  and  4th  order  reduced 
model  at  small  voltage  input  V  =  O.luoKs.  The  center  point  displacement  is  plotted  against  time. 

With  the  mesh  size  of  40  x  20,  an  880th  order  system  was  generated  by  the  linearization  process,  i.e.  880 
ODE’s.  We  compare  the  time  history  of  the  displacement  of  the  center  point  of  the  beam  at  different  step 
input  voltages:  0.1V,  2V  and  9V.  Three  models  are  compared,  namely,  the  full  finite  difference  solution  of 
the  original  non-linear  model  of  equations,  the  linearized  system  of  equations,  and  a  4th  order  reduced  model 
generated  using  a  Krylov-subspace  method  for  selecting  Vq.  Figure  8  shows  that  for  a  small  input  voltage 
V  -  0.1  volts  the  three  curves  representing  solutions  of  each  of  the  three  models  overlap  with  one  another. 
The  figure  shows  that  with  such  a  small  input  voltage  the  original  system  behaves  almost  perfectly  linearly, 
and  that  the  4th  order  reduced  model  faithfully  reproduces  the  behavior  of  the  880th  order  linear  system. 

In  Figure  9,  we  see  that  the  linearized  model  starts  to  deviate  from  the  non-linear  model,  but  the  4th  order 
reduced  model  still  follows  the  linear  model  nearly  exactly.  Figure  10  demonstrates  that  at  the  pull-in  voltage, 
the  time  response  of  the  structure  is  extremely  non-linear.  The  linear  model  and  the  4th  order  macro-model 
are  only  accurate  during  the  initial  part  of  the  transient. 

For  a  discretization  size  of  10  x  5,  which  is  in  turn  a  70th  order  linear  system,  we  compared  the  frequency 
responses  of  the  linear  model  and  various  orders  of  Arnoldi  based  macro-models.  Figure  11  shows  the 
comparison  of  the  frequency  responses  of  the  large  linear  system  and  two  macro-models  that  are  of  the  order 
of  2  and  10,  respectively.  We  see  from  Figure  11  that  the  original  linear  system  is  a  well  damped  system.  The 
original  linear  system  has  a  bandwidth  frequency  of  1.8  x  105. 

This  figure  also  shows  that  the  2nd  order  macro-model  perfectly  matches  the  linear  model  in  a  low  frequency 
range  up  to  10®  H z.  And  the  10th  order  model  is  able  to  follow  all  the  oscillatory  behavior  both  in  the  gain 
plot  and  the  phase  plot.  The  frequency-domain  accuracy  of  the  10th  order  model  would  be  important  if  the 
device  where  part  of  a  feedback  system. 

5  A  Circuit  Representation  for  Micromachined  devices 

The  semi-analytical  macromodeling  approach  can  be  quite  effective  for  design  synthesis  and  optimization,  but 
given  the  macromodel  has  built-in  assumptions  about  device  behavior,  the  approach  is  not  a  very  effective 
verification  strategy.  The  numerical  model-order  reduction  approach,  even  if  the  difficulties  with  nonlinearities 
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V=2  volts 


Figure  9:  At  V  =  2 volts,  the  linear  model  starts  to  deviate  from  the  non-linear  model.  The  reduced-order 
model  still  follows  the  linear  model. 


V=9  volts,  non-linear  pull-in  occurs 


Figure  10:  A  comparison  of  the  time  responses  of  the  non-linear  model,  linear  model  and  4th  order  reduced 
model  at  small  voltage  input  V  =  0.1  volts.  The  center  point  deflection  is  plotted  against  time. 


Frequency  Response  Comparison 


Figure  11:  A  comparison  of  the  frequency  responses  of  the  full  linear  system,  2nd  order  reduced  model  and 
10th  order  reduced  model. 


were  overcome,  seems  poorly  suited  to  the  synthesis  and  optimization  phase  of  a  design  because  the  approach 
requires  a  complete  layout  of  the  device,  and  it  yields  no  information  about  sensitivities  to  changes  in  design 
attributes.  In  this  section  we  will  describe  a  third  alternative  to  extending  circuit  simulation  to  include 
micromachined  devices.  An  approach  will  be  developed  which  comes  much  closer  to  providing  a  schematic-like 
description  for  micromachined  devices,  but  at  the  cost  of  narrowing  the  range  of  micromachined  devices  which 
can  be  so  treated.  So  in  that  sense,  this  circuit-like  description  intended  to  make  simulation  easier  is  also  a 
step  towards  top-down  or  structured  design  methodologies  [46,  47,  48,  49,  50,  51,  52,  53,  54], 

Developing  a  circuit  representation  for  micromachined  devices  involves  determining  the  list  of  elements 
for  the  circuit  representation,  the  model  for  each  element,  and  the  definition  of  the  nature  or  discipline 
for  the  terminals  of  the  elements.  The  goals  in  selecting  the  list,  model  and  nature  include  design  reuse  and 
simulation  accuracy.  Element  parameterization  provides  both,  while  supporting  a  wide  class  of  micromachined 
device  designs.  Parameterization  with  both  design  attributes  and  process  parameters  (captured  in  the  model 
technology  file)  allows  process  independent  models  that  can  be  used  to  simulate  devices  in  a  variety  of 
fabrication  lines.  A  conservative  Kirchhoffian  network  representation  is  used  both  for  simulation  accuracy,  and 
for  compatibility  with  electronics  design.  Signal-flow  representations,  commonly  used  for  behavioral  or  system- 
level  modeling,  are  more  cumbersome  to  use  for  this  application  because  they  are  based  on  unidirectional 
elements  while  mechanical  and  electrical  components  interact  bidirectionally. 

5.1  Element  Hierarchy 

Micromachining  technology  combines  sacrificial  etching  with  VLSI-style  deposit,  pattern  and  etch  sequences 
to  produce  miniaturized  mechanical  components  that  are  suspended,  cavitied,  hinged,  or  otherwise  mounted. 
The  circuit  approach  described  below  focuses  primarily  on  suspended  microstructural  devices  as  that  is  the 
most  mature  micromachining  design  space.  In  principle,  the  circuit  approach  can  be  extended  to  include 
hinged  structures  for  optical  applications  or  cavitied  structures  for  fluidic  applications. 

Suspended  micromachined  devices  involve  plates  tethered  by  beams  to  anchors.  Air  gaps  between 
conductive  micromachined  elements  act  as  a  variable  sensing  capacitance  and  a  source  for  electrostatic 
actuation  force.  For  example,  plates  can  be  actuated  electrostatically  to  tilt  micromirrors  in  digital  computer 
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Figure  12:  Atomic  elements  for  design  of  suspended  micromachined  systems  include  (a)  anchor,  (b)  beam,  (c) 
gap  and  (d)  plate. 


Figure  13:  (a)  Layout  and  (b)  Schematic  Representation  of  an  individual  crab-leg  microresonator;  and  (c) 
Layout  and  (d)  Schematic  representation  of  “O”  coupling  spring. 


projection  displays  [2]  or  in  optical  switches.  Micromachined  inertial  sensors  employ  one  or  more  plates  as 
proof  masses,  which  move  when  accelerated  and  whose  motion  is  sensed  capacitively  [55].  Such  suspended 
micromachined  structures  decompose  into  anchors,  beams,  plates  and  gaps,  as  shown  parametrized  in 
Figure  12.  This  set  of  elements  is  chosen  for  three  reasons:  they  all  occur  commonly  (albeit  sized  by  appropriate 
geometric  parameters);  they  are  modular  (in  the  sense  that  they  are  decoupled  from  neighboring  elements); 
and  their  behavior  can  be  accurately  approximated  with  a  simple  lumped  parameter  model.  For  the  class  of 
Manhattan-geometry  suspended  polysilicon  devices,  this  set  of  four  elements  completely  covers  the  possible 
design  space,  and  forms  a  set  of  atomic  basis  elements  in  the  circuit  representation. 

A  circuit  simulation  environment  for  micromachined  devices  based  on  this  element  library  with 
parametrized  behavioral  models,  called  NODAS  (Nodal  Design  of  Actuators  and  Sensors),  has  been  devel¬ 
oped  [56,  57,  60].  Schematic  examples  of  a  crab-leg  resonator  and  an  O-shaped  spring  using  the  basis 
elements  are  shown  in  Figure  13.  The  use  of  circuit  element  libraries  for  nodal  analysis  of  micromachined 
systems  [64,  65,  66]  and  for  microgyroscope  simulation  [58,  67]  are  also  being  simultaneously  pursued. 

Circuit  representation  of  additional  elements  at  higher  levels  of  the  hierarchy  may  also  be  desirable, 
primarily  because  such  elements  aid  in  the  capture  of  complex  designs.  In  particular,  the  parameterized 
functional  elements  such  as  the  linear  comb-drive  sensor  or  actuator,  or  the  crab-leg  spring,  O-spring,  or 
folded-flexure  spring  are  easily  re-used  because  they  capture  a  single  function  (generate  electrostatic  force, 
provide  mechanical  stiffness,  etc.)  and  hence  can  be  accurately  represented  by  behavioral  models.  While 
parametrized  models  at  high  levels  of  component  abstraction  are  still  possible,  the  fixed  topology  of  these 
components  limits  their  re-usability.  Moreover,  the  large  number  of  design  variables  in  such  abstractions 
significantly  increases  the  complexity  of  generating  a  parametrized  lumped- parameter  model  as  detailed  in 
Section  3.2. 

In  addition  to  carefully  choosing  the  elements  in  the  design  library  to  ensure  richness  of  coverage  of  the 
design  space,  the  terminals  of  each  element  needs  to  be  carefully  chosen  to  balance  the  need  for  interoperability 
between  the  elements  and  the  accuracy  and  speed  of  design  simulation.  By  using  the  same  terminal  natures 
at  all  levels  of  the  design  hierarchy,  a  composable  design  representation  for  mixed-level  simulation  is  possible. 
This  is  particularly  important  as  simulation  of  entire  systems  at  the  atomic  level,  though  possible,  may  require 
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unnecessary  long  simulation  times.  A  library  consisting  of  the  most  common  atomic  and  functional  elements 
therefore  supports  both  rapid  simulation  as  well  as  the  capability  to  represent  a  wide  class  of  designs. 

As  the  underlying  simulation  representation  is  a  Kirchhoffian  network,  the  nature  of  the  quantity  defined 
across  and  through  each  branch  in  the  network  is  very  important.  In  the  electrical  domain,  voltage  across 
and  current  through  a  branch  is  the  accepted  standard.  For  mechanical  domains,  no  standard  nature  exists. 
Two  possible  translational  mechanical  across-through  relations  are  velocity-force  and  displacement-force.  The 
latter  representation  is  preferred  in  micromechanica!  design,  as  displacement  is  generally  the  most  common 
observable  state.  Similarly,  the  rotational  mechanical  nature  is  angular  displacement  across  and  torque  through 
a  branch. 

The  associated  reference  directions  of  the  mechanical  terminals  correspond  directly  to  the  physical 
directions  of  displacement  and  force.  As  with  electrical  circuit  simulation,  a  consistent  and  systematic  set  of 
associated  reference  directions  for  mechanical  terminals  is  essential.  A  simple  convention  specifies  translational 
displacement  across  variables  as  positive  in  the  positive-axis  directions  and  the  rotational  displacement 
variables  as  positive  in  a  counterclockwise  rotation  (right-hand  rule)  around  the  positive-axis  directions. 
Through  variables  going  into  a  branch  are  interpreted  as  positive- valued  force  or  torque  in  the  positive  direction 
[60]. 

Once  the  terminal  natures  are  defined,  the  element  can  be  modeled  by  relating  the  flow  through  the 
terminals  to  the  potential  across  the  terminals.  This  model  is  often  called  a  constitutive  relationship  in 
network  theory.  The  models  need  to  capture  all  the  physics  of  the  given  element,  hence  a  beam  element  needs 
to  include  mass,  spring  and  damping  physics,  all  parametrized  by  the  beam  design  geometry  and  the  process 
model  parameters.  Parametrized  models  that  are  within  a  few  percent  of  continuum  simulations  have  been 
derived  using  techniques  described  in  Section  3. 

Mechanical  parasitics  need  to  be  considered  for  accurate  circuit  simulation.  For  example,  due  to  the  lumped 
parameter  modeling  of  the  atomic  elements,  the  joint  between  two  beams  in  a  flexure  becomes  a  parasitic. 
The  compliance  of  the  joint  is  a  fringing  effect  that  can  be  modeled  by  extending  the  length  of  the  beams 
incident  at  the  joint.  If  one  of  the  beams  incident  at  the  joint  is  significantly  wider  than  the  other  beam, 
then  the  moment  relations  at  the  joint  need  to  be  considered.  Extension  factors  and  the  use  of  plate  joints 
have  been  verified  by  comparing  the  circuit  simulation  with  continuum  finite  element  simulation  for  all  the 
common  flexure  topologies  and  a  range  of  beam  sizes.  In  all  cases,  the  error  in  flexure  compliance  and  resonant 
frequency  was  less  than  2%.  Additional  sources  of  parasitics  include  the  capacitance  and  resistances  in  the 
interface  between  the  microstructures  and  electronics. 

5.2  Micromechanical  Bandpass  Filter  Circuit 

The  flexibility  of  the  microelelectromechanical  circuit  representation  is  best  demonstrated  by  returning  to  the 
bandpass  filter  example.  The  filter  is  composed  of  three  identical  resonators,  each  with  a  center  frequency  of 
300  kHz,  coupled  by  beam  springs  [9].  The  topology  of  the  filter  (Figure  14)  with  both  mechanical  structures 
and  interface  circuitry  is  captured  in  the  schematic  using  the  symbols  from  the  NODAS  element  library  [68, 61]. 
The  interface  circuitry  includes  Q-adjustment,  frequency  tuning,  and  a  trans-resistance  sense  amplifier. 

When  an  a.c.  input  voltage,  Vin,  Ls  applied  across  the  electrostatic  comb  drive,  the  suspended  shuttle 
masses  and  flexural  beams  will  be  driven  by  the  electrostatic  force  and  move  in  the  x  direction.  This  mechanical 
vibration  is  coupled  to  the  other  two  resonators  via  the  coupling  beams,  resulting  in  three  resonant  peaks, 
thus  forming  a  pass  band.  The  location  and  spacing  of  the  three  peaks  are  determined  by  the  stiffness  of 
coupling  beams,  leading  to  different  center  frequency  and  bandwidth.  An  equivalent  SPICE  representation  is 
derived  in  [9,  68]  using  the  methods  described  in  Section  3.2.  The  equivalent  SPICE  models  represent  the 
mechanical  resonators  as  second-order  systems  of  lumped  parameters  for  mass,  spring  constant  and  damping, 
and  represent  the  coupling  beams  as  massless  ideal  springs  with  a  coupling  spring  constant.  Figure  15  shows 
the  result  of  simulating  the  filter  in  NODAS  with  massless  beam  models  compared  to  the  equivalent  SPICE 
model  in  vacuum.  The  natural  frequency  of  the  resonators  is  299.43 kHz,  and  the  quality  factor  is  495,000. 
The  coupling  beams  are  88.2 g,m  long  and  1.12 fim  wide.  There  are  three  peaks  around  the  natural  frequency, 
ranging  from  299 AZkHz  to  299.95 kHz.  NODAS  and  SPICE  results  match  to  within  4%. 

The  peaks  can  be  flattened  to  form  a  flat  passband  by  applying  Q-adjustment  series  resistors,  shown  by 
the  simulated  filter  frequency  response  in  Figure  16.  The  three  sharp  peaks  of  the  initial  high-Q  filter  are  now 
compressed  down  to  a  nearly  flat  passband  with  a  ripple  of  -21  dB  (Q  of  587).  NODAS  and  SPICE  simulation 
results  match  to  within  4%. 

The  actual  flexure  and  coupling  beams  have  finite  mass,  and  can  not  be  treated  as  ideal  springs,  a 
simplification  needed  for  the  analytical  derivation  of  the  equivalent  SPICE  model.  Flexure  beam  mass  shifts 
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Figure  14:  Circuit  Schematic  of  three  resonator  filter 


the  center  frequency  of  the  filter,  while  coupling  beam  mass  contributes  to  the  lumped  parameter  equivalent 
masses  of  adjacent  resonators  also  shifts  resonant  frequencies  as  well  as  causes  passband  distortion.  This 
combined  effect  can  be  quantified  by  comparing  the  frequency  response  using  the  default  NODAS  beam  model 
(with  finite  mass)  with  the  response  using  a  massless  beam  model,  as  shown  in  Figure  17. 

The  combination  of  ease  of  schematic  entry,  simulation  accuracy,  applicability  in  iterative  design  of  the 
wide  class  of  suspended  microsystems,  compatibility  with  existing  VLSI  design  flows,  and  support  for  co¬ 
simulation  of  electronic  and  micromachined  devices  make  this  micromecahnical  circuit  simulation  approach 
very  attractive.  To  expand  this  circuit  approach,  continued  research  is  needed  in  identifying  the  basis  elements 
for  enlarged  design  spaces  that  include  cavitied  and  hinged  structures.  Terminal  natures  for  additional  physics 
such  as  fluidic  pressure  and  flow  rate  or  optical  beam  intensity  are  also  needed  to  model  new  classes  of  devices. 
Additionally,  methodologies  and  tools  for  automated  extraction  of  geometric  and  material  parameters  for 
accurate  simulation  are  crucial  for  the  wide  applicability  of  this  simulation-based  design  approach. 

5.3  Extraction  from  Layout 

The  circuit-like  representation  for  micromachined  devices  fits  perfectly  with  the  synthesize  and  optimize 
phase  of  device  design,  as  device  performance  can  be  simulated  but  layout  details  can  be  avoided.  For 
this  representation  to  be  useful  during  the  verification  phase,  it  must  be  possible,  to  extract  the  circuit  from 
the  layout.  Just  like  for  mainstream  integrated  circuits,  layout  extraction  involves  recognizing  patterns  in  the 
layout  and  then  inferring  a  one-to-one  correspondence  between  the  layout  patterns  and  the  circuit  elements. 

Layout  extraction  involves  recognition  of  the  layout  patterns  that  correspond  to  the  circuit  schematic 
elements  based  on  their  features  (shape,  size,  location).  Once  the  schematic  elements  are  recognized,  the 
extraction  creates  a  connected  schematic  to  capture  the  shape  and  location,  and  annotates  the  element  sizes, 
thus  creating  a  complete  schematic  representation  of  the  microstructure  layout.  The  elements  can  be  extracted 
as  fixed  values  (e.g.,  plate  has  lug  mass),  or  as  geometrical  parameters  (e.g.,  square  plate  has  length  of  100/i.m). 
The  abstractions  used  for  mixed-domain  circuit  simulation  are  based  on  geometrically  parametrized  models  of 
the  atomic  elements,  requiring  extraction  of  geometrical  parameters  from  the  layout.  This  approach  is  similar 
to  device  extraction  in  VLSI,  where  geometrical  parameters  for  the  MOS  model  is  extracted  from  the  layout. 
Unlike  VLSI  layout  extraction,  however,  the  features  (shape,  size  and  position)  of  each  layout  rectangle  is  of 
utmost  importance  in  recognizing  the  constitutive  micromachined  atomic  elements  (VLSI  extractors  would 
consider  a  sequence  of  beams  forming  a  suspention  to  be  a  single  wire). 

Once  the  constitutive  atomic  elements  are  recognized,  element-specific  extraction  can  be  used  as  necessary. 
This  procedure  involves  purely  geometrical  reasoning  to  identify  the  constitutive  schematic  elements,  followed 
by  determining  the  appropriate  parameters  for  each  instance  of  an  atomic  or  functional  element  [63]  found 
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Figure  18:  Folded  flexure  resonator  (a)  layout,  (b)  canonical  representation,  (c)  canonical  representation  after 
separating  the  fingers,  (d)  intermediate  state,  (e)  detected  state,  (f)  functional  elements  detected 


in  the  layout.  In  addition  to  extracting  a  schematic  representation  from  the  layout,  parasitics  can  also  be 
identified  and  extracted  [53,  54]. 

Rectangles  that  comprise  the  layout  are  generated  by  algorithms  specific  to  the  layout  editing  tools. 
The  first  step  in  feature  recognition  for  Manhattan-geometry  structures  therefore  involves  creating  a  unique 
representation  of  the  layout.  Starting  from  an  input  layout  in  CIF  (Caltech  Interchange  Form,  shown  in 
Figure  18(a)),  the  rectangles  in  the  layout  axe  partitioned  into  a  canonical  representation,  such  that  each 
rectangle  (or  cell)  has  only  one  neighbor  on  each  side  (shown  in  Figure  18(b)).  The  use  of  the  canonical 
representation  allows  the  development  of  algorithms  that  are  independent  of  the  CAD  software  used  to  generate 
the  input  layout.  The  disadvantage  of  the  canonical  representation  is  a  significant  increase  in  the  number  of 
rectangles  to  be  processed.  Most  of  this  increase  comes  from  the  presence  of  fingers  in  the  design,  hence 
they  are  removed,  and  the  layout  re-canonized,  as  shown  in  Figure  18(c).  The  functionality  of  each  of  the 
cells  is  then  determined  by  its  shape,  size  and  connectivity.  Non-structural  mask  layers  (such  as  those  that 
define  anchors)  are  used  to  obtain  hints  for  possible  functional  uses  for  each  of  the  cells,  using  rules  from 
a  process  description  file.  Also  contained  in  this  file  are  rules  for  atomic  element  recognition,  for  example, 
cells  with  one  connected  side  are  cantilever  beam  fingers,  and  cells  with  connections  on  opposing  sides  are 
considered  to  be  beams.  The  partitioning  due  to  the  canonical  representation  algorithm  results  in  multiple 
adjacent  cells  performing  the  same  function.  These  multiple  cells  have  to  be  combined  to  minimize  the  number 
of  unnecessary  nodes  in  the  netlist.  Cell  merging,  first  in  the  horizontal  direction,  and  then  in  the  vertical 
direction  accomplishes  this  for  the  mass  and  anchor  cells.  This  merging  reduces  the  total  number  of  ports  in 
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the  generated  netlist,  hence  contributes  to  the  management  of  the  simulation  time  for  the  extracted  netlist. 
The  resulting  netlist  directly  corresponds  to  the  atomic  elements  (beams,  plates,  gaps  and  anchors)  as  shown 
in  Figure  18(e)  in  the  circuit  representation  of  Section  5.1.  Thus  the  primary  objective  of  having  a  check  on 
the  designed  layout  can  be  achieved.  Device  function  can  also  be  confirmed  via  the  “circuit”  simulation  in 
Section  5.2  on  this  netlist. 

Higher-level  functional  element  models  can  be  detected  by  processing  the  extracted  netlist.  A  functional 
element  library  containing  rules  for  detecting  commonly  used  spring  suspensions  and  comb-drive  topologies 
is  external  to  the  extraction  tool,  and  can  be  customized  to  alternate  processes  and  design  styles.  Finger 
orientation,  region  of  occurrence,  and  geometrical  parameters  (length,  width  and  inter-finger  gap)  are  used 
to  partition  the  set  of  recognized  fingers,  which  are  then  analyzed  for  connectivity  resulting  in  the  extracted 
comb-drives.  Spring  detection  is  accomplished  via  a  finite  state  machine  (FSM)  based  algorithm.  Starting  from 
a  start  state  (always  an  anchor  atomic  element),  the  type  of  beam  and  joint  determines  transitions  into  the 
intermediate  states,  and  onto  the  final  state,  which  indicates  the  type  of  spring  detected.  The  joint  transitions 
are  classified  according  to  the  number  of  ports  and  the  direction  of  rotation,  and  provide  the  fundamental 
abstraction  on  which  this  FSM-based  detection  works.  The  FSM  for  each  of  the  springs  is  described  in  the 
component  library.  The  connected  sets  of  beams  and  springs  obtained  after  the  atomic  recognition  is  passed 
through  each  of  these  FSMs  to  recognize  their  type,  as  shown  by  the  example  in  Figure  18(f).  Simulation- 
based  verification  using  this  level  of  extraction  is  an  order  of  magnitude  faster  than  at  the  atomic  element 
level,  and  is  seen  to  be  crucial  for  an  iterative  design  methodology. 

The  challenge  for  the  extraction  methodology  is  to  provide  a  rich  set  of  basic  recognition  functions  and 
a  language  for  combining  these  functions  in  both  process-independent  element  recognition,  and  process- 
dependent  use  of  layer  information  to  support  recognition  and  extraction  of  all  layouts  that  can  be  mapped  to 
the  circuit  element  library.  This  has  yet  to  be  demonstrated  for  non-Manhattan  geometry  structures  in  single- 
structural  layer  polysilicon  micromachining  processes.  However,  for  the  commonly  used  Manhattan-geometry 
suspended  microstructure  design  style  and  polysilicon  micromachining  process,  extraction  has  been  extremely 
effective  in  microstructure  layout  verification. 


6  Conclusions 

In  this  survey  paper  we  presented  and  contrasted  three  different  approaches  for  extending  circuit  simulation 
to  include  micromachined  devices.  The  most  commonly  used  method,  that  of  using  physical  insight  to  develop 
parameterized  macromodels,  is  presented  first.  The  issues  associated  with  fitting  the  parameters  to  simulation 
data  while  incorporating  design  attribute  dependencies  was  shown  to  require  sophisticated  intervention.  In 
addition,  the  semi-analytic  approach  did  not  seem  to  provide  a  very  effective  verification  path.  Then,  the 
numerical  model  order  reduction  approach  to  macromodeling  was  presented,  and  it  was  shown  that  the  key 
difficulty  remains  finding  automatic  way  to  perform  nonlinear  model  reduction.  In  addition,  model-order 
reduction  seemed  to  be  ineffective  during  the  synthesis  and  optimization  phase  of  design,  because  no  attribute 
sensitivities  were  computed.  Lastly,  we  described  the  recently  developed  circuit-based  approach  for  simulating 
micromachined  devices,  and  described  the  design  hierarchy  and  the  use  of  a  catalog  of  parts.  We  also  showed 
that  the  circuit-based  approach  can  be  combined  with  extraction  from  layout,  providing  an  effective  approach 
to  verification.  The  only  short-coming  of  this  circuit-based  approach  is  that  only  some  design  styles  and 
technology  can  be  supported. 
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